#include "cubicRoot.hpp"
#include "typeAliases.hpp"

#include <catch2/catch_test_macros.hpp>
#include <catch2/catch_approx.hpp>

#include <boost/math/tools/cubic_roots.hpp>
#include <iostream>

using Catch::Approx;
using namespace hms;
using namespace hms::cubic;
using boost::math::tools::cubic_roots;

TEST_CASE("Depressed cubic coeffs"){
	SECTION("General case, fundamental data types"){
		scalar
			a { 2},
			b {-6},
			c {-1},
			d { 4};
		REQUIRE( cubic::p(a,b,c  ) == Approx(-3.5) );
		REQUIRE( cubic::q(a,b,c,d) == Approx(-0.5) );
	}
	SECTION("Simplified case, fundamental data types"){
		scalar
			a { 1},
			b {-6},
			c { 0},
			d { 4};
		REQUIRE( cubic::p(a,b,c  ) == Approx(-12) );
		REQUIRE( cubic::q(a,b,c,d) == Approx(-12) );
		REQUIRE( cubic::norm_noLin::p(b  ) == Approx(-12) );
		REQUIRE( cubic::norm_noLin::q(b,d) == Approx(-12) );
	}
	SECTION("General case, Eigen types"){
		Array2s
			a { 2,-1},
			b {-6, 3},
			c {-1, 2},
			d { 4,-3};
		REQUIRE( cubic::p(a,b,c  ).isApprox( Array2s{-3.5, -5} ) );
		REQUIRE( cubic::q(a,b,c,d).isApprox( Array2s{-0.5, -1} ) );
	}
	SECTION("Simplified case, Eigen types"){
		Array2s
			a { 1, 1},
			b {-6, 3},
			c { 0, 0},
			d { 4,-3};
		REQUIRE( cubic::p(a,b,c  ).isApprox( Array2s{-12, -3} ) );
		REQUIRE( cubic::q(a,b,c,d).isApprox( Array2s{-12, -1} ) );
		REQUIRE( cubic::norm_noLin::p(b  ).isApprox( Array2s{-12, -3} ) );
		REQUIRE( cubic::norm_noLin::q(b,d).isApprox( Array2s{-12, -1} ) );
	}
}

TEST_CASE("Discriminant"){
	SECTION("Fundamental data types"){
		scalar
			a { 2},
			b {-6},
			c { 2},
			d { 4},
			p, q;
		auto set_pq = [&](){
			p = cubic::p(a,b,c  );
			q = cubic::q(a,b,c,d);
		};
		set_pq();
		REQUIRE( cubic::discriminant(p,q) == Approx(-5) );
		b = -4;
		set_pq();
		REQUIRE( cubic::discriminant(p,q) == Approx(116) );
	}
	SECTION("Eigen data types"){
		Array2s
			a { 2, 2},
			b {-6,-4},
			c { 2, 2},
			d { 4, 4},
			p { cubic::p(a,b,c  ) },
			q { cubic::q(a,b,c,d) };
		REQUIRE( cubic::discriminant(p,q).isApprox( Array2s{-5, 116} ) );
	}
}

TEST_CASE("cubic root computation, casus irreducibilis"){
	auto setGeneral = [](
		const auto& a, const auto& b, const auto& c, const auto& d,
		auto& p, auto& q, auto& scale, auto& acosTerm,
		auto& t0, auto& t1, auto& t2,
		auto& r0, auto& r1, auto& r2
	){
		p = cubic::p(a,b,c  );
		q = cubic::q(a,b,c,d);
		scale = threeReal::scale(p);
		acosTerm = threeReal::acosTerm(p,q);
		t0 = threeReal::depressedRoot(scale, acosTerm, 0);
		t1 = threeReal::depressedRoot(scale, acosTerm, 1);
		t2 = threeReal::depressedRoot(scale, acosTerm, 2);
		r0 = cubic::root(t0, a, b);
		r1 = cubic::root(t1, a, b);
		r2 = cubic::root(t2, a, b);
	};
	auto setSimplified = [](
		const auto& b, const auto& d,
		auto& p, auto& q, auto& scale, auto& acosTerm,
		auto& t0, auto& t1, auto& t2,
		auto& r0, auto& r1, auto& r2
	){
		p = cubic::norm_noLin::p(b  );
		q = cubic::norm_noLin::q(b,d);
		scale = threeReal::scale(p);
		acosTerm = threeReal::acosTerm(p,q);
		t0 = threeReal::depressedRoot(scale, acosTerm, 0);
		t1 = threeReal::depressedRoot(scale, acosTerm, 1);
		t2 = threeReal::depressedRoot(scale, acosTerm, 2);
		r0 = cubic::norm_noLin::root(t0, b);
		r1 = cubic::norm_noLin::root(t1, b);
		r2 = cubic::norm_noLin::root(t2, b);
	};
	SECTION("Fundamental data types"){
		scalar
			a { 2},
			b {-6},
			c { 2},
			d { 4},
			p, q,
			scale, acosTerm,
			t0, t1, t2,
			r0, r1, r2;
		SECTION("General case"){
			setGeneral(a,b,c,d,p,q,scale,acosTerm,t0,t1,t2,r0,r1,r2);
			REQUIRE( cubic::discriminant(p,q) < 0 );
			REQUIRE( scale == Approx(1.633) );
			REQUIRE( acosTerm == Approx(0.91174) );
			scalar res { 0.618034 };
			REQUIRE( t0 == Approx(1) );
			REQUIRE( t1 == Approx( res) );
			REQUIRE( t2 == Approx(-(1+res) ) );
			REQUIRE( r0 == Approx(2) );
			REQUIRE( r1 == Approx( 1+res) );
			REQUIRE( r2 == Approx(-res) );
			/* compare against boost implementation */
			std::array<scalar,3> br { cubic_roots(a,b,c,d) };
			// std::cout
			// 	<<   "hms   roots: " << r0 << ", " << r1 << ", " << r2
			// 	<< "\nboost roots: " << br[0] << ", " << br[1] << ", " << br[2]
			// 	<< '\n';
			REQUIRE( r0 == Approx(br[2]) );
			REQUIRE( r1 == Approx(br[1]) );
			REQUIRE( r2 == Approx(br[0]) );
		}
		SECTION("Simplified case"){
			setSimplified(b,d, p,q, scale,acosTerm, t0,t1,t2, r0,r1,r2);
			REQUIRE( cubic::discriminant(p,q) < 0 );
			REQUIRE( scale == Approx(4) );
			REQUIRE( acosTerm == Approx(0.240911416) );
			scalar
				res0 {  0.88448370193933  },
				res1 {  0.884250603336952 },
				res2 { -0.768734305276285 };
			REQUIRE( t0 == Approx( res0+3) );
			REQUIRE( t1 == Approx( res1-2) );
			REQUIRE( t2 == Approx( res2-2) );
			REQUIRE( r0 == Approx( res0+5) );
			REQUIRE( r1 == Approx( res1  ) );
			REQUIRE( r2 == Approx( res2  ) );
			/* compare against boost implementation */
			std::array<scalar,3> br { cubic_roots(1.,b,0.,d) };
			REQUIRE( r0 == Approx(br[2]) );
			REQUIRE( r1 == Approx(br[1]) );
			REQUIRE( r2 == Approx(br[0]) );
		}
	}
	SECTION("Eigen data types"){
		Array2s
			a { 2,-1},
			b {-6, 3},
			c {-1, 2},
			d { 4,-3},
			p, q,
			scale, acosTerm,
			t0, t1, t2,
			r0, r1, r2;
		SECTION("General case"){
			setGeneral(a,b,c,d,p,q,scale,acosTerm,t0,t1,t2,r0,r1,r2);
			REQUIRE( (cubic::discriminant(p,q) < 0).all() );
			REQUIRE( scale   .isApprox( Array2s{2.16024689946929, 2.58198889747161} ) );
			REQUIRE( acosTerm.isApprox( Array2s{0.45702710293778, 0.44542446588745} ) );
			
			REQUIRE( t0.isApprox( Array2s{ 1.93853719123054, 2.33005873956798} ) );
			REQUIRE( t1.isApprox( Array2s{-0.14370504950751,-0.20163967572340} ) );
			REQUIRE( t2.isApprox( Array2s{-1.79483214172303,-2.12841906384458} ) );
			REQUIRE( r0.isApprox( Array2s{ 2.93853719123054, 3.33005873956798} ) );
			REQUIRE( r1.isApprox( Array2s{ 0.85629495049249, 0.79836032427660} ) );
			REQUIRE( r2.isApprox( Array2s{-0.79483214172303,-1.12841906384458} ) );
		}
		SECTION("Simplified case"){
			setSimplified(b,d, p,q, scale,acosTerm, t0,t1,t2, r0,r1,r2);
			REQUIRE( (cubic::discriminant(p,q) < 0).all() );
			REQUIRE( scale   .isApprox( Array2s{4, 2} ) );
			REQUIRE( acosTerm.isApprox( Array2s{0.240911415937805, 0.349065850398866} ) );
			Array2s
				res0 {  0.88448370193933 ,  0.87938524157182 },
				res1 {  0.884250603336952, -0.34729635533386 },
				res2 { -0.768734305276285, -0.53208888623796 };
			REQUIRE( t0.isApprox( res0 + Array2s{ 3, 1} ) );
			REQUIRE( t1.isApprox( res1 + Array2s{-2, 0} ) );
			REQUIRE( t2.isApprox( res2 + Array2s{-2,-1} ) );
			REQUIRE( r0.isApprox( res0 + Array2s{ 5, 0} ) );
			REQUIRE( r1.isApprox( res1 + Array2s{ 0,-1} ) );
			REQUIRE( r2.isApprox( res2 + Array2s{ 0,-2} ) );
		}
	}
}

TEST_CASE("Cubic root computation, single real root"){
	using Eigen::sqrt;
	using std::sqrt;
	SECTION("Fundamental data types"){
		scalar p, q, sqrtTerm, c1, c2, t;
		p = -2;
		q = -4;
		REQUIRE( cubic::discriminant(p,q) > 0 );
		sqrtTerm = cubic::singleReal::sqrtTerm(p,q);
		c1 = cubic::singleReal::cubeRoot<0>(q, sqrtTerm);
		c2 = cubic::singleReal::cubeRoot<1>(q, sqrtTerm);
		t  = cubic::singleReal::depressedRoot(q, sqrtTerm);

		REQUIRE( sqrtTerm == Approx(sqrt(3.) * 10 / 9) );
		REQUIRE( c1 == Approx(1.577350269) );
		REQUIRE( c2 == Approx(0.4226497308) );
		REQUIRE( t  == Approx(2) );
	}
	SECTION("Eigen data types"){
		Array2s
			p {-2, -3},
			q {-4, -6},
			sqrtTerm, c1, c2, t;
		REQUIRE( (cubic::discriminant(p,q) > 0).all() );
		sqrtTerm = cubic::singleReal::sqrtTerm(p,q);
		c1 = cubic::singleReal::cubeRoot<0>(q, sqrtTerm);
		c2 = cubic::singleReal::cubeRoot<1>(q, sqrtTerm);
		t  = cubic::singleReal::depressedRoot(q, sqrtTerm);

		REQUIRE( sqrtTerm.isApprox(Array2s{sqrt(3.) * 10 / 9, 2*sqrt(2)}) );
		REQUIRE( c1.isApprox(Array2s{1.57735026918963 , 1.799632345152}) );
		REQUIRE( c2.isApprox(Array2s{0.422649730810374, 0.555669052456122}) );
		REQUIRE( t .isApprox(Array2s{2, 2.35530139760812}) );
	}
	SECTION("Normalised, no linear term"){
		Array2s
			b {-sqrt(6),-3},
			d {-4-(-2*pow(6.,1.5))/27,-4},
			r;
		ArrayNNs<2,cubic::norm_noLin::bufCols> buf;
		r = cubic::norm_noLin::root(b,d,buf);
		REQUIRE( r.isApprox(Array2s{2.81649658092773, 3.35530139760812}) );
		/* compare against boost implementation */
		REQUIRE( r[0] == Approx(cubic_roots(1.,b[0],0.,d[0])[0]) );
		REQUIRE( r[1] == Approx(cubic_roots(1.,b[1],0.,d[1])[0]) );
	}
}